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We present a newly enhanced version of the Monte Carlo Shell Model method by in- 
corporating the conjugate gradient method and energy-variance extrapolation. This new 
method enables us to perform large-scale shell-model calculations that the direct diagonal- 
ization method cannot reach. This new generation framework of the MCSM provides us 
with a powerful tool to perform most-advanced large-scale shell-model calculations on cur- 
rent massively parallel computers such as the K computer. We discuss the validity of this 
O ' method in ab initio calculations of light nuclei, and propose a new method to describe the 

^ , intrinsic wave function in terms of the shell-model picture. We also apply this new MCSM to 

■ the study of neutron-rich Cr and Ni isotopes using the conventional shell-model calculations 
with an inert *°Ca core and discuss how the magicity of A*' = 28, 40, 50 remains or is broken. 

§1. Introduction 

^: 

\ The understanding of the structure of ah nuclei from the first principle, called 

t:;^- ■ usually the ab initio nuclear calculation, is one of the ultimate goals in nuclear 

theory. For this purpose, one usually starts with the nucleon degree of freedom, 
\ i.e., protons and neutrons as the building block of a nucleus, assuming the free 

■ (or bare) nucleon-nucleon force to be the interaction between nucleons. Recent ah 
initio calculations often include not only the two-nucleon force but also the three- 
nucleon force which is specific to the interaction between composite particles like 

^ ■ nucleons. Several ab initio methods have been proposed, and their accuracy has 

^ , been investigated in great detail, for instance, in terms of the binding energy of the 

four-nucleon system.^^ It remains, however, rather difficult or infeasible to go beyond 
systems with A > 12, where A is the number of nucleons. This is largely due to 
strong nucleon-nucleon correlations that requires a large number of single-particle 
states to be included for the description of many-body states in ground states or 
low excited states. Although this problem can be resolved to a certain extent by 
using effective interactions based on various renormalization techniques, the number 
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of many-body states to be included remains prohibitively large in most cases, and 
increases exponentially with the nucleon number. In this paper, we present a new 
version of the Monte Carlo shell model (MCSM), demonstrating how it can contribute 
to ab initio nuclear calculation as well as to conventional but quite huge shell model 
calculations. 

The nuclear shell model has been known to be successful in describing the struc- 
ture of atomic nuclei based on nuclear forces. It dates back to 1949 when Mayer and 
Jensen discovered the magic numbers as a consequence of a mean potential including 
the spin-orbit coupling.'^' While the original concept of the shell model of Mayer and 
Jensen is basically an independent particle picture, the concept has been modified 
and extended significantly over decades since then. The modern shell model uses a 
sufficiently large number of many-body basis states which are superposed utilizing 
properly constructed effective interactions so as to provide us with accurate many- 
body eigenstates. The many-body basis state is a Slater determinant usually with 
single-particle states taken from a harmonic oscillator potential. 

In conventional shell-model calculations, an inert core is assumed: all single- 
particle states below a given magic number are completely occupied, forming a closed 
shell. Single-particle states between this magic number and the next magic number 
constitute a valence shell, and nucleons in this shell are called valence nucleons. In 
the conventional shell model, only valence nucleons are activated. The single-particle 
states of activated nucleons in the shell-model calculation are called the model space. 
The model space is a concept for calculation, and is the same as some valence space 
in many cases. But, in other cases, the model space can be taken wider or smaller 
than the relevant valence shell depending on some interest or limitation. We thus 
distinguish model space from valence shell hereafter. Note that the model space is 
a more computational concept. 

An effective interaction is obtained for valence nucleons, and is defined for each 
model space. Effects from the inert core, those from virtual excitations from the inert 
core and those from virtual excitations to states above the model space are assumed 
to all be incorporated into this effective interaction by renormalizing it appropriately. 
This can be a very important issue, but the conventional shell model assumes that 
there is such an interaction, while its determination can be phenomenological. 

It has been shown that many nuclear properties with 8 < N[Z) < 20 are excel- 
lently described with the conventional shell model calculation in the Is-Od (model) 
space (often called the sd shell) ,'3® and those with 20 < A^(Z) < 32 are also quite 
well described in the Ip-Of (model) space (often called the pf shell) .1^® Note that 
the effective interactions used in those calculations are hybrid products of micro- 
scopic derivation and empirical adjustments. 

While the conventional shell model assumes a relatively limited model space as 
exemplified above, the same computational procedure is applicable to the ab initio 
calculation when no inert core is assumed and a large number of single-particle states 
are taken so that the calculation becomes close to calculations in the whole Hilbert 
space. Here, we refer to this ab initio method within the shell model as the ab initio 
shell model, one of which is the no-core shell model (NCSM),'^ a famous model with 
great success. 
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Whether the ab initio shell model or the conventional shell model (with a core) 
is used, all one has to do in computation is to diagonalize a Hamiltonian matrix 
spanned by all the possible many-body states in a given model space (i.e., single- 
particle space of activated nucleons). The number of relevant many-body basis states 
determines the dimension of this matrix. This dimension is often called the shell- 
model dimension, and causes a serious computational issue as we shall see later. The 
number of the many-body states is NpCup x NnCn„ without symmetry consideration, 
where Np {N^) is the number of proton (neutron) single-particle states taken and Up 
(un) is the number of protons (neutrons) activated. It roughly increases exponen- 
tially with Np (or A^^^) and also with Up (or n„). Hence, without even considering 
the ab initio shell model with a huge number of Np and Nn, the conventional shell 
model for heavier nuclei already suffers from the huge dimensionality of the Hamil- 
tonian matrix necessary to tackle the whole nuclear chart because the number of 
valence orbits {Np or Nn) increases for heavier nuclei. For instance, the dimension 
needed for A ^ 80 N = Z nuclei is estimated to be ~ 10^^ without any symmetry 
consideratiorP' when the lp-0f-2s-ld-0g orbits are assumed to be the valence shell. 
Although this number can be reduced by 1-2 orders of magnitude by choosing only 
the states of the same Jz as taken usually (i.e., M-scheme calculation), it is still far 
from the current computational limit of 10^*^^^^ in the M-scheme. 

In order to go beyond the computational limit of the shell model, a new method 
for performing the shell-model calculation named the Monte Carlo shell model (MCSM) 
has been developed since 1996,^' guided by the Quantum Monte Carlo Diagonal- 
ization method.'^ The MCSM utilizes the idea of the auxiliary field Monte Carlo 
method that is taken in the Shell Model Monte Carlo method,!^ but this method and 
the MCSM are completely different. The MCSM aims to represent many-body states 
with a small number of highly selected many-body basis states. In this sense, the 
MCSM is regarded as an "importance truncation" of the entire many-body space.'^ 
The basis state should be represented in a compact form (i.e., a wave function with 
a small number of parameters), it should be able to approximate the nuclear many- 
body state efficiently, and its matrix elements should be calculated easily. To meet 
this demand, we usually use deformed Slater determinants, and parity and total 
angular momentum projection are operated on them if needed. Otherwise, a pair 
condensed stat^^^ or a quasi-particle vacuum state, which is used in the Hartree- 
Fock-Bogoliubov calculation and also taken as the basis state of the VAMPIR (abbr. 
for Variational After Mean field Projection In Realistic model spaces) method,^ is 
a good candidate when required. Note that the VAMPIR method has been devel- 
oped with more emphasis on calculating the energy spectra rather than improving 
the energy of a specified state. The MCSM basis states are added one by one, and 
each basis state to be added is selected among many candidate Slater determinants 
generated stochastically so that the energy of the state under consideration can be as 
low as possible. This step is repeated until the energy converges sufficiently. Thus, 
this original MCSM is characterized as "stochastic" in terms of the method of basis 
variation, and as "sequential" in terms of the procedure of basis variation. 

From the computational point of view, the MCSM has an advantage over the 
conventional diagonalization method in tolerance to the increase of the model space 
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and the particles. When it is assumed that the number of basis states needed in 
a MCSM wave function is fixed, the total computational cost is scaled by the cost 
of each Hamiltonian matrix element. This is roughly proportional to {Nsp)" with 
a ~ 3-4 in the case of Nsp = Np = Nn, being much milder than the exponential 
increase. This advantage had enabled us to perform the full p/-shell calculation in 
^^Ni (with ~ 10^ M-scheme dimension) with good accuracj^^ several years before 
the exact diagonalization was carried out.'^ See the review papei'^ for more details 
of early achievements. 

Recently, the computational environment has been changing rapidly. The num- 
ber of available CPU cores is expanding to be typically in the range of tens of 
thousands for the world leading supercomputers as compared to a few tens to hun- 
dreds of CPU cores used in the early MCSM calculations. The K computer will 
contain more than 700,000 cores upon its completion in the autumn of 2012. This 
situation has strongly motivated us to renew the MCSM method to be suitable using 
up-to-date massive parallel computers, in addition to wider applications including 
ab initio shell model calculations. Since 2009, we have developed a new MCSM 
package, which includes not only renewed code, but also a novel methodology and 
numerical algorithm.^iS Among the methodological advancements, the most im- 
portant is the introduction of the energy-variance extrapolation method.^ In the 
original MCSM, the energy of a many-body state is evaluated directly from the en- 
ergy expectation value of the MCSM wave function, which must be higher than the 
exact solution. It was quite hard to know the difference between this value and the 
exact value. The energy-variance extrapolation method serves as a powerful tool 
to pin down accurately the exact solution from a series of approximated solutions. 
Another methodological change is incorporating a variational aspect into the MCSM 
by varying the basis state, which enhances the lowering of calculated energy values. 
Equipped with the variational-type improvement, the MCSM is now characterized as 
"deterministic" in terms of the method of basis variation. Thus, keeping the original 
idea, the present MCSM, associated with several advancements, can be regarded as 
a new generation. 

In this paper, we review the outline of the new-generation MCSM and show 
some of its earliest applications that were not feasible with the original MCSM. This 
paper is organized as follows. In Sect. [21 the outline of the new-generation MCSM 
and its feasibility are presented. We also discuss the adaptability of the MCSM 
to massively parallel supercomputers. In Sect. [21 application to the ab initio shell 
model is demonstrated. Along with numerical success, the intrinsic shape and the 
clustering of light nuclei are also discussed with the MCSM wave function. In Sect. 
[H application to the neutron-rich chromium and nickel isotopes is presented as a case 
of medium-heavy mass nuclei. This region is being intensively studied in radioactive 
isotope facilities over the world, and also challenges nuclear models because several 
shapes coexist and are mixed, which results in the rapid change of the dominant 
shape over the isotopes. In Sect. \5\ we summarize this paper and give an outlook 
and future perspective of the MCSM toward the launch of the K computer. 

We note here that many new features may appear in the structure of exotic 
nuclei, because unbalanced numbers of protons and neutrons may create situations 
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where unknown or little known aspects of nuclear forces become visible and pro- 
duce large impacts. The shell evolution due to the tensor (oic^^'^^'^^ and the 
three-body forc^^D are some examples. Note that the shell evolution explains basic 
trends of single-particle properties and one needs comprehensive calculations to ob- 
tain physical quantities and look into correlations in depth. The exploration of such 
unknown features need theoretical framework directly linked to nuclear forces, and 
we expect a large contribution from the new generation of MCSM. 

§2. Outline of the new-generation Monte Carlo shell model 

The new-generation MCSM can be divided roughly into two stages: the varia- 
tional procedure to obtain the approximated wave function and the energy-variance 
extrapolation utilizing the obtained wave function. We briefly describe these two 
parts in Sects. 12.21 and 12. 3|, respectively. The shell-model Hamiltonian and the form 
of the variational wave function used in the MCSM is shown in Sect 12. II The ad- 
ditional improvement to make the extrapolation procedure stable is demonstrated 
in Sect. 12.41 The numerical aspects of the framework mainly for massively parallel 
computation are discussed in Sect. 12. 5[ 

2.1. Shell-model Hamiltonian and variational wave function 

In conventional nuclear shell-model calculations assuming an inert core, we use 
a general two-body interaction: 

H = H^^^ + H^^^ = J2ticUi+ Yl ^ijkicWjCiCk, (2-1) 

i i<j,k<l 

where c| denotes a creation operator of single particle state i. is a one-body 

Hamiltonian with single-particle energies tj, and iJ^^^ is a two-body interaction which 
has parity and rotational symmetry and is represented by so-called Two-Body Matrix 
Elements (TBMEs).l2a 

In the case of ab initio shell-model calculations, the Hamiltonian is taken as 

H = T - TcM + ^ = %4ci + X] Vijkic\c]ciCk, (2-2) 

ij i<j,k<l 

where T and Tcm are the total kinetic energy and the kinetic energy of the center-of- 
mass motion, respectively. Note that Tcm has both one-body and two-body compo- 
nents. The V represents two-nucleon interaction, e.g., the JISP16 interactiorP^ in 
Sect. 3. In this work, because we do not treat explicitly three-nucleon forces, both 
the Hamiltonians of these two kinds of shell-model calculations consist of one-body 
and two-body interactions. 

If necessary, the removal of spurious center-of-mass motion can be done by uti- 
lizing the prescription of Gloeckner and LawsoEpSJ to suppress the contamination of 
the center-of-mass motion. In this prescription the Hamiltonian to be diagonalized 
is replaced by 

H' = H + /3cm-f^cm (2-3) 
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with Hera being defined as 

Hem = + -MAco^Yi^ - -hw, (2-4) 

2AM 2 2 ' ^ ^ 

where R and P are the position and momentum of the center of mass. By taking 
/3cm large enough, (Hem) is suppressed as a small value. 

In the framework of the MCSM, the approximated wave function is written as a 
linear combination of angular-momentum- and parity-projected Slater determinants, 

Nb I 

i'^iv,) = x; E &pmM, (2-5) 

n=l K=~I 

where A^;, is the number of the Slater-determinant basis states. The Plfj^ operator 
is the angular- momentum and parity projector defined as 

,7. _ l + 7r77 2/ + l r ^ r 



PUk = ^ J dQ DW{Q)e^'^'^e^^'ye^'^'% (2-6) 

where 17 = (a, /3, 7) are the Euler angles and D^j^p^{Q) denotes Wigner's D-function. 
n stands for the parity transformation. Each |(/)„) is a deformed Slater determinant 
defined in Eq. (jA-ip . The coefficients f^^^ are determined by the diagonalization 
of the Hamiltonian matrix in the subspace spanned by the projected Slater deter- 
minants, PMx\4>n) with —I < K < I and 1 < n < A*";,. This diagonalization also 
determines the energy, E^^ = {'F]yjH\'Fiy^), as a function of A^^. Note that the 
dimension of the subspace is (21 -|- 1)A'';,, not Nj,. The Slater determinant basis, 
is given by variational calculation to minimize EN^=n- We increase N), until 
E^^ converges enough, or the extrapolated energy converges. The strategy of this 
variational calculation will be discussed in the next subsection. 

2.2. Variational procedure 

In this section, we discuss the efficient process for the determination of Z)(") in 
Eq. (jA-ip to minimize E]sf^ and the history of its developments in this section. In the 
original MCSM calculation, the basis states are selected from many candidates gen- 
erated stochastically utilizing the auxiliary-field Monte Carlo technique, the detail 
of which was discussed in Ref . fT2|) . In order to determine in Eq.dZS]) efficiently, 
M. Honma et al. introduced the few-dimensional basis approximation,^®' in which 
they adopted the steepest gradient method. It succeeded in estimating the energy 
of p/-shell nuclei with a relatively small number of bases, however, the direction 
of the gradient is not necessarily the most efficient choice to reach the local mini- 
mum and it is practically difficult to control the step width in the gradient direction. 
Mainly because of this problem, the number of basis states was limited to be rather 
small (iVfo ~ 30) compared to that in the MCSM (iVt ~ 100). On the other hand, 
W. Schmid et al. adopted the quasi-Newton method for energy minimization in 
the VAMPIR approximation to obtain the optimized quasi-particle vacua as basis 
states.^^' G. Puddu demonstrated that the quasi-Newton method works also in the 
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Slater determinant bases. In the quasi-Newton method, the step width is automati- 
cally determined to minimize the energy along the modified gradient directionP^'EHli 
As a result, the step width is relatively large in early steps of the quasi-Newton pro- 
cedure along the sophisticated gradient direction so that the additional basis state 
has sufficiently linearly independent component of the other basis states. 

In the present work, we adopt the conjugate gradient (CG) method,'^''^ which 
also includes the linear minimization in the modified gradient direction. In compari- 
son to the quasi-Newton method, the CG method is expected to be advantageous to 
save memory usage and improve computational efficiency for parallel computation 
because a large Hessian matrix is utilized in the quasi-Newton method, which is not 
used in the CG method. Note that the energy of the linear combination of projected 
Slater determinants is optimized, not that of the unprojected Slater determinants. In 
this sense, this variational procedure is "variation after projection and configuration 
mixing" . 

We proposed the following four steps to optimize the parameters of deformed 
Slater determinants: 

1. Monte Carlo sampling utilizing auxiliary field technique (original MCSM), 

2. Sequential optimization for each basis state (SCG), 

3. Refinement process of each basis states repeatedly (refinement), 

4. Full (simultaneous) optimization of all basis states (FCC). 

In the first step, we perform the original MCSM procedure to obtain the approx- 
imated wave function, which can be used as an initial state of the CG process. The 
detail of the original MCSM is not discussed here, but in Ref. In Refs. fT2]) .[T7 11 . 
we select the best basis states from an order of 1000 candidates, which is generated 
stochastically. However, the necessary number of candidates can be far suppressed 
if we proceed to the next step. 

In the second step, called the SCG method, we perform a variational calculation, 
with variational parameters being D^^\ sequentially by minimizing the Ei without 
changing the other bases, 

Z)(2),...d(»-i), which were already fixed. In other 
words, we perform variational calculations with a set of variational parameters, 
sequentially. We increase the number of basis states, A'^^, until the energy reaches 
sufficient convergence. 

In the third step which is called the refinement process, we take an initial state 
from the result of the SCG calculation. In the SCG calculation, the D^^') is not the 
best optimized parameters to minimize E^^, since -D^^^ is determined to minimize 
El. Then, we first fix the number of basis states, Nm, and restart to minimize the 
energy Ejy^ by the CG method to optimize each basis state from D^^^ to one 
by one. We iterate a few times to perform a routine of the CG process for all basis 
states one by one to get a better energy calculation. 

In the fourth step called the FCG, we fix at the beginning, and determine 
aU coefficients D^^b) ^ith 1 < Afo < iV„ by minimizing the E^^ at once. In other 
words, we perform variational calculation with a set of all variational parameters, 
Jji^b) ^ simultaneously to minimize the E^^. Its initial state is generated by the 
refinement process in order to save the amount of computation time. In principle, 
the FCG yields the closest energy to the exact one within a fixed number of basis 
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states, Nra, and the refinement process also reaches the energy provided by the FCG 
with a large number of iterations. The refinement process needs a far smaller amount 
of time than that in the FCG, and provides us with an approximation good enough 
to the solution in the FCG. 

Figure [1] shows the convergence pattern of the ground-state energy of ^^Ge in 
/sPffg shell, which consists of the O/7/2, 1^3/2; lPi/2) and Og'9/2 orbits. Its M-scheme 
dimension is very large and amounts to 140,050,484. However, it can be handled 
by the recent shell-model diagonalization code.^^-^ The energy eigenvalue E^^^ = 
{^n^\H\1'nC) is plotted against the number of basis states, Ni,. The SCG method 
attains faster convergence than the original MCSM: the SCG gives the same energy as 
the original MCSM in almost half the number of basis dimension. The SCG method 
enables us to compute the energy variance in a smaller amount of the computation 
time than that of the original MCSM, since the time to compute the energy variance 
is proportional to Nm{Nm + l)/2. 

In the case of the FCG, the number 
of basis states of the variational wave 
function is fixed to Nm = 32 and we plot 
energy expectation values in the sub- 
space spanned by the N^, basis states. In 
the FCG aU basis states, D^^\D^'^\ 
and D(32) 

are optimized to minimize 
£^32, while in the SCG wave function 
D^^^ is optimized to minimize Ei. Thus, 
the calculated Ei of the FCG is much 
higher than the Ei of the SCG and 
the El of the original MCSM. However, 
£^32 of the FCG is lower than that of 
the SCG and the original MCSM, which 
means the FCG provides us with the 
best approximation for a fixed number 
of basis states. Therefore, we should discuss the energy convergence as a function of 
Nm in the FCG, not as a function of A'^;,, unlike the case of the SCG. 

While the FCG yields the best energy with the fixed Nm, the FCG iteration 
needs much computational resources since we have to replace all basis states at 
every FCG iteration. Which level of the optimization is most efficient in terms of 
the computation time depends on the property of the objective wave function. For 
applications, we adopt the original MCSM method and the SCG method in Sect. [3] 
and the refinement process in Sect. |H 

2.3. Energy-variance extrapolation 

Even if the variational procedure works excellently, a small gap between the 
exact energy eigenvalue and the energy expectation value of the approximated wave 
function remains. The gap can be removed by extrapolation procedures, which have 
been studied intensively.'SS'Kli'ESli'ESli Among them, the energy-variance extrapola- 
tion method is a general framework for the supplementation of the variational cal- 
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Fig. 1. Convergence patterns of the ground- 
state energy of '^^Ge in the fspgg shell. The 
circles, open triangles, and squares denote 
the resuhs of the original MCSM, SCG, 
and FCG with Nm = 32, respectively. 
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culation and is rather independent of the form of the approximated wave function. 
This method has been known in condensed matter physics-—' and was firstly intro- 
duced in shell-model calculations with conventional particle-hole truncation by T. 
Mizusaki and M. Imada!^ In the present work, we apply this extrapolation method 
to estimate the exact eigenvalue precisely by utilizing a sequence of the linear com- 
binations of the projected Slater determinants, I^^at^) with 1 < A'^^ < which 
have been obtained in the variational procedure discussed in Sect l2.2[ The major 
obstacle of the energy-variance extrapolation is the necessity of the computation of 
the energy variance. Its efficient computation is described in Sect l2.5l 

The energy-variance extrapolation method is based on the fact that the en- 
ergy variance of the exact eigenstate is zero. The energy variance of the approx- 
imated wave function is not exactly zero, but rather small and approaches zero 
as the approximation is improved. In the framework of the energy-variance ex- 
trapolation, we draw the energy i^TVj = {'^Nt\H\'I^Nt) against the energy variance 
{AH^)Ni = i^NjH^l'l'Nt) - E% , the plot of which is called "energy- variance plot" . 
The variance usually approaches zero as Ni^ increases, as the point in the energy- 
variance plot approaches the y-intercept. Following the idea of Ref j34i) . These points 
are usually fitted by a first- or second-order polynomial such as 

E = co + ci{AH^) + C2{{AH^)f, (2-7) 

where these coefficients cq, ci, and C2 are determined by least square fit. By ex- 
trapolating the fitted curve into the y-intercept we obtain the extrapolated energy, 
namely, cq. 

In the framework of the present study, the variational procedure discussed in 
Sect. 12.21 provides us with a sequence of approximated wave functions, which can be 
utilized in the energy-variance extrapolation method. The SCG procedure provides 
us with a successive sequence of the wave functions, I^n^) with 1 < A^fe < Nm, 
where Nm is the maximum of Nf,. For each Nh, we evaluate the energy E]sf^ and 
energy variance {AH'^)^,,- Here, we demonstrate how the extrapolation method 
works with ^^Ni in the pf shell as an example. The effective interaction is the FPD6 
interaction,'^^ which was adopted also in Ref. [T7]) . 

Figure [2)[a) shows the energy-variance plot of the SCG method for yrast states 
of ^^Ni. In this calculation, we take the K = I state only in the angular- momentum 
projector, namely /^^^ = if K ^ I in Eq. (|2-5p . for simplicity. As Nf, increases the 
energy-variance point moves smoothly and approaches the y-axis or variance zero, 
except for the S"*" state. The fitted curves for these points are shown as red solid lines, 
and they also show smooth behavior. The extrapolated energies, or y-intercepts of 
the fitted curves, agree quite well with the exact ones, which are shown as open 
circles on the y-axis. Especially concerning the ground-state energy, the minimum 
variance of the approximated wave function is {AH'^) ^,,=100 = 0.89MeV^, which is 
smaller and gives a better approximation than the result of the original MCSM, 
{AH'^)n^=150 = imUeV'^^ mainly thanks to the introduction of the CO method. 

On the other hand, the plot of the S"*" state shows the anomalous behavior, 
in which the energy decreases when increasing the number of basis states but the 
variance does not decrease at {AH'^) ~ 4MeV^. As a result, the simple extrapolation 



10 




10 
Energy variance (MeV^) 



J=10 



J=8 



J=6 



-5=4 



J=0 



20 

Energy variance (IVIeV^) 



40 



Fig. 2. Energy-variance plot for ^^Ni (a) without reordering, (b) with reordering. The open circles 
on the y-axis denote the exact shell-model energies. See text for details. 



method apparently fails. A straightforward solution to this failure is to increase Nf, 
until the extrapolation method works, and it is shown in Ref. ^2S). However, in 
other cases, it might be difficult due to the increase of computation time. We discuss 
another remedy for such cases in the following section. 

2.4. Reordering technique to improve the energy-variance extrapolation 

The anomalous behavior can be removed by the reordering of the basis states.'^ 
In this section, we demonstrate that the reordering technique makes the energy- 
variance extrapolation stable and avoids the difficulty of an anomalous kink such as 
the 8^ state of ^^Ni discussed in the previous subsection. 

A sequence of the approximated wave functions, {^2), \'^Nt)j •••) I'^A'm)) 

is specified by a set of basis states and its order (I \(j)2), I'pN,,), I'pNm))- In the 
SCG method, the order of the basis states is determined by the variational procedure. 
However, we can shuffle the order of basis states and make another sequence of the 
approximated wave functions without additional computational effort. If we assume 
that the extrapolated value is independent of the order of the basis states, there exists 
the best order which makes the extrapolation procedure stable. In the reordering 
method, the order is determined so that the gradient of the curve in the energy- 
variance plot is as small as possible. The practical algorithm to determine the order 
is discussed in Ref. [39]) . 

Figure [2{b) shows the extrapolation 
plot with the reordering technique using 
the same set of the basis states as in Fig. 
^a). In this case, the gradient of the fit- 
ted curve is so stable that the points in ^ 
the energy-variance plot are fitted by a 2 
first-order polynomial and the region to i "i^i 

UJ 

-182 

exact 
J=0*i < 

'"""^^O 10 20 30 

Energy Variance (MeV ) 

Fig. 3. Energy vs. energy variance plot of the 
ground state energy of ^'^Ge in the fspgg 
shell. The wave functions are provided 
by the FCG with A'',^ = 8 (triangles), 16 
(squares), 24 (diamonds), and 32 (circles). 
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be used for the fit is rather broad. In 
Fig. [2]^b), the anomalous behavior of 
the variance plot of the 8"*" state disap- 
pears and the extrapolated value agrees 
with the exact energy quite well. ^^Ni is 
also known to have a shape coexistence 
featurej^S which is plausible to be the 
origin of the anomalous behavior of the 
8"^ state in Fig. [2|^a) like the case of 
^^Ge discussed in Refl39]). 

We demonstrate another example 
of the energy-variance extrapolation 

method combined with the reordering technique in Fig. [3] using the FCG wave 
functions of ^^Ge with the JUN45 effective inter action.HD' In this figure, the en- 
ergy expectation value of the FCG wave function is plotted as a function of the 
corresponding energy variance with the basis states being Nm = 8, 16, 24, and 32, 
and their second-order fitted curves. The order of the basis states in each sequence 
is determined by the reordering technique.'^S Note that there is no specific order 
in the FCG procedure because we treat all basis states on an equal footing. The 
extrapolated energy apparently converges except for Nm = 8. While we show the 
second-order fit in the figure, even the extrapolated value of the first-order fit agrees 
well with that of the second-order fit in the case of Nm = 32. 

2.5. Computational aspects 

The large-scale shell-model calculation is one of the challenging issues for nu- 
clear physics. It is essential to develop a program which runs efficiently on recent 
supercomputers. Concerning a calculation using a single CPU, we proposed a sophis- 
ticated way of efficient computation of the matrix element of non-orthogonal Slater 
determinants in Ref l42p . Moreover, since the main trend of recent supercomputers 
favors massively parallel computers, the parallel efficiency is worth discussing for 
future studies. 

In the case of the SCG method, we need to compute the Hamiltonian matrix 
elements and the gradient vector concerning D*^") of two angular- momentum-, parity- 
projected Slater determinants. The three-dimensional integral of the Euler angles in 
Eq. (|2-6p is performed by discretizing each range of the integral into mesh points using 
the Gauss-Chebyshev quadrature for the z-axis rotations and the Gauss-Legendre 
quadrature for the y-axis rotation, which is shown in Eq. ()A-3p . The numbers of the 
mesh points are taken as 26 for the z-axis rotation and 16 for the y-axis rotation for 
example. The parity projection is equivalent to two mesh points. The product of 
these mesh points, which are denoted by A in Eq. ()A-3p give rise to A^mesh = 26^ x 
16 X 2 = 21632 components in evaluating a matrix element in terms of the projected 
Slater determinants. Since these components can be computed independently, the 
program was written for massive parallel computation. When we apply the matrix- 
product technique discussed in RefH2|) with the bunch size Abunch being e.g. 30, we 
still have the A^[,-^mesh/-^bunch — 721A^fe elements to be computed in parallel. 
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Figured] shows the parallel efficiency of the benchmark calculation of the ground 
state of ^^Ge as an example. The model space consists of the pf shell and gg/2 
orbit and the PFG9B3 effective interaction is used:^^' Its M-scheme dimension 
reaches 1.7 x 10^^, which is far beyond the current limitation of the Lanczos method. 
The MCSM result of this system was already reported in Refs. [T7j) . [39j) . This 
benchmark was performed using the Intel Fortran compiler ver. 11. cP on the T2K 
open Supercomputer at the University of Tokyo.l^ 




Number of cores Number of cores 



Fig. 4. (a) Speedup of the parallel computation of the SCG process in unit of the computation time 
using 16 CPU cores. The squares, open circles, triangles and filled circles represent the inverses 
of computation times of the variational process to obtain 1st, 4th, 16th, and 32nd basis states 
respectively. The solid line shows ideal scaling to guide the eyes, (b) Speedup of the parallel 
computation of the energy variance. The squares, triangles, and circles represent the inverses of 
computation times of 1, 4, and 8 basis states, respectively. 

Figure S] (a) shows the performance gain of the parallel computation of the 
SCG process to determine the 1st, 4th, 8th, 16th, and 32nd basis states respectively. 
Although the parallel efficiency for calculating the first basis state is not good because 
of the small amount of computation, the efficiency for the 32nd basis state with 2048 
CPU cores reaches 82% of that with 16 cores. 

We calculate the energy variance using the formula shown in AppxEP^ Be- 
cause the two-body matrix elements in the M-scheme, Vijki, are sparse due to the 
symmetry which the Hamiltonian has, we store in memory only non-zero matrix 
elements in block-diagonal form by treating f as a rank-2 matrix with indices 

{ij) and (kl). Thus, we can compute the energy variance efficiently, and the detail 
of the practical computation is written in Appx|Xl In a similar manner to the case 
of the variational process, we compute the energy variance by dividing the whole 
computation into matrix elements which are moreover divided into each mesh point 
of Eq. ()2-6p . resulting in Nh{Ni, + 1) x A'mesh independent components to be computed 
in parallel. In addition, we do not need an iterative process like the CG method, and 
therefore a small amount of network communication appears only at the beginning 
and at the end of the computation. Thus the performance scaling of the parallel 
computation seems perfect at Nf, > 4, which is shown in Fig. IHb). 

In practice, it took totally 807 seconds to obtain an SCG wave function of ^^Ge 
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0^ state with 32 basis states, and it took 588 seconds to compute the energy variance 
of this SCG wave function using 2048 CPU cores. 

§3. Application of the MCSM to the ab initio shell model 

In this section we focus on the latest appUcation of the MCSM to the ab initio 
shell model calculations, which has become feasible recently with the aid of major 
development of the MCSM algorithm discussed in Sect. [2] and also a remarkable 
growth in the computational power of state-of-the-art supercomputers. First, the 
no-core shell model (NCSM) and its variants are briefly reviewed. The limitation of 
the NCSM and the motivation for the application of the MCSM to the ab initio no- 
core Full Configuration Interaction (FCI) approach are further discussed here. Then, 
the current status of the benchmarks in the no-core MCSM is referred to based 
on the results mostly from Ref.SS' Finally our challenge to visualize the intrinsic 
states constructed by superpositioned non-orthogonal Slater determinants is also 
demonstrated. 

3.1. Ab initio shell models 

One of the major challenges in nuclear theory is to understand nuclear structure 
and reactions from ab initio methods. Ab initio calculations for nuclear many-body 
systems beyond A = 4 have recently become feasible due to the rapid evolution 
of computational technologies these days. In ab initio approaches for the nuclear 
structure calculations, all the nucleons constituting the nucleus are considered as 
the fundamental degrees of freedom and the bare/effective interactions based on 
realistic nuclear forces are adopted. As for bare two- and three-nucleon interac- 
tions, the phase-shift equivalent family of two-nucleon interactions, derived from 
the meson-exchange theory and chiral Effective Field Theory, in addition to three- 
nucleon interactiond23'SH)'llSli.[SSli.[5D jg usually used. 

Ab initio NCSM has been emerging for about a decade and is now available 
for the study of nuclear structure and reactions in the p-shell nuclei.^ Unlike the 
conventional shell model with a core, the NCSM does not assume an inert core just 
like the name itself implies and treats all the nucleons composing the nucleus on 
an equal footing. The NCSM is thus said to be one of the ab initio approaches 
along with the Green's Function Monte CarloS2J a^d Coupled Cluster theory.lSB In 
the NCSM (in a narrow sense^ the model (or basis) space is usually truncated by 
the so-called Nmax, which is the sum of the excitation quanta above the reference 
state. The effective interactions renormalized to that model space are used so as to 
obtain the faster convergences of the energy with respect to Nmax- Generally, the 
effective interactions are derived by the so-called Lee-Suzuki-Okamoto method.'^S 
The NCSM result approaches the exact solution either by taking the larger model 
space with the level of the cluster approximation fixed or by improving the order of 
the cluster expansion with the model space fixed. 

A similar but distinct approach to the NCSM is the No-Core Full Configuration 
(NCFC) approach.ESJ The NCSM result by using the effective interactions derived by 
the Lee-Suzuki-Okamoto procedure approaches the exact solution either from below 
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or above due to the violation of the strict variational upper bound of the exact solu- 
tion. Therefore the extrapolation of the NCSM result into the infinite model space is 
obscure. The NCFC method employs the bare or effective (softened) low-momentum 
interactions evolved from bare nuclear forces by the renormalization group transfor- 
mations!^ which validates the variational upper bound of the calculated energy. 
The NCFC enables access to full ab initio solutions by a simple extrapolation into 
the infinite model space in the two-dimensional parameter space {Hlo, Nmax)- One of 
the advantages both in the NCSM and NCFC methods is the perfect factorization of 
the intrinsic and the center-of-mass wave functions, so that the intrinsic state does 
not suffer from the spurious center-of-mass motion. 

As ah initio approaches treat all of the nucleons democratically, computational 
demands for the calculations explode exponentially as the number of nucleons and/or 
the model spaces increase. Current limitation of the direct diagonalization of the 
Hamiltonian matrix by the Lanczos iteration is around the order of 10^'^ shown in 
Fig. [5j So far the largest calculations have been done in the ^"^N with N^ax = 8 
which results in the M-scheme many-body matrix dimensions being ~ 10^ and as- 
sociated non- vanishing three- nucleon force matrix elements being ~ 4 x 10^^.^^' In 
order to access heavier nuclei beyond the p-shell region with larger model spaces 
by ab initio shell-model methods, many efforts have been devoted for several years. 
One of these approaches in the Nmax truncation is the Importance- Truncated NCSM 
(IT-NCSM).IS3 In the IT-NCSM, the model spaces are extended by using the im- 
portance measure evaluated by the perturbation theory. Another approach is the 
Symmetry- Adapted NCSM (SA-NCSM),ES where the model spaces are truncated 
by the selected symmetry groups. 

Besides the Nmax truncation of the 
model space in the ah initio shell mod- 
els, there is the FCI method to give 
the exact solutions in the fixed model 
space. Different from the N^ax trun- 
cation in the NCSM and NCFC meth- 
ods, the FCI truncates the model space 
by the single-particle states, so-called 
N shell or emax{= N shell " !)• As shown 
in Fig. [5l the explosion of the dimen- 
sionality prohibits the full ah initio so- 
lutions of the FCI (and also the NCSM) 
beyond the lower p-shell region. Sim- 
ilar to the attempts of the IT-NCSM 
and SA-NCSM, the MCSM is one of the 
promising candidates to go beyond the 
FCI method.'^''^ Note that there is a similar approach to the no-core MCSM re- 
ferred to as the Hybrid Multi-Determinant method.'® In the following subsection 
we will show some recent investigations by the ah initio no-core MCSM. 
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3.2. Benchmarks of the MCSM to the ah initio no-core FCI 

As an exploratory work of the original MCSM has been applied to the no- 
core calculations for the structure and spectroscopy of the beryllium isotopes.^ In 
Ref.l^ the low-lying excited states of ^'^Be and ^^Be are investigated. The excitation 
energies of the first and second 2"*" states and the B(E2; 2^ — )■ 0^^ ) for ^"Be with a 
treatment of spurious center-of-mass motion show good agreement with experimental 
data. The deformation properties of the 2^ and 2^ states for ^"Be and of the 2^ state 
for ^^Be are studied in terms of electric quadrupole moments, E2 transitions and the 
single-particle occupations. The triaxial deformation of ^'^Be is also discussed in 
terms of the B(E2; 2^ — )• 2^) value. This work motivates a further extension of the 
MCSM application to the ah initio FCI calculations .'^SJ Currently, the availability of 
the MCSM for the no-core calculations has been tested extensively in light nuclei.'^ 
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Fig. 6. *He ground-state energies as functions of number of basis states (left) and energy variance 
(right). The red, green, blue and purple solid symbols (horizontal dashed hues in the left figure 
and open symbols at the zero energy variance in the right figure) are the MCSM (FCI) results 
in Ngheii = 2, 3, 4 and 5, respectively. The harmonic oscillator energies are taken at optimal 
values for each state and model space. The Coulomb interaction and the spurious center-of-mass 
motion effect are not considered. Isospin symmetry is assumed. 

As a typical example, the behavior of the ground-state energies of ^He (0"*") with 
respect to the number of basis states and to the energy variance in N shell = 2 — 5 
are shown in Fig. [6l Figure [7| illustrates the comparisons of the energies for each 
state and model space between the MCSM and FCI methods. The FCI gives the 
exact energies in the fixed size of the mode space, while the MCSM gives approxi- 
mated ones. Thus the comparisons between them show how well the MCSM works 
in no-core calculations. For this benchmark comparison, the JISP16 two-nucleon 
interaction is adopted and the Coulomb force is turned off. Isospin symmetry is 
assumed. The energies are evaluated for the optimal harmonic oscillator frequencies 
where the calculated energies are minimized for each state and model space. Here 
the contributions from the spurious center-of-mass motion are ignored for simplic- 
ity. In Fig. [71 the comparisons are made for the states; ^IIe(0"^), ^He(0"'"), ^Li(l"'"), 
^Li(l/2-, 3/2-), SBe(0+), i°B(l+, 3+) and i2c(0+). The model space ranges from 
N shell = 2 through 5 for A < 6 (4 for A>1). Note that the energies of i°B(l+, 3+) 
and ^^C(O^) in N shell = 4 ^ne available only from the MCSM results. The M-scheme 
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Fig. 7. Comparisons of the energies between the MCSM and FCI along with the fully converged 
NCFC results where availableP^ The MCSM (FCI) results are shown as the solid (dashed) lines 
that nearly coincide where both are available. The extrapolated MCSM results are illustrated by 
bands. From top to bottom, the truncation of the model space is Nshm = 2 (red), 3 (green), 4 
(blue) and 5 (purple). Note that the MCSM results are extrapolated by the energy variance with 
the second-order polynomials.'^ Also note that all of the results of and ^^C at Ngheii = 4 
were obtained only with MCSM. 

dimensions for these states (1.82 x 10^° for M = I and 1.52 x 10^° for M = 3 in 
^"^B and 5.87 x 10^^ for M = in ^^C) are already marginal or exceed the current 
limitation in the FCI approach. The number of basis states are taken up to 100 in 
^sheii = 2 — 4 and 50 in Ngheii = 5. In Fig. [71 the solid (dashed) hues indicate the 
MCSM (FCI) results. The shaded regions express the extrapolations in the MCSM, 
and the lower bound of the shaded region corresponds to the extrapolated energy. 
Furthermore, we also plot the NCFC results for the states of 4 < A < 8 as the fully 
converged energies in the infinite model space. As seen in Fig. [71 the energies are 
consistent with each other where the FCI results are available to within ~ 100 keV 
(~ 500 keV) at most of the MCSM results with(out) the energy-variance extrapola- 
tion in the MCSM. The other observables besides the energies also give reasonable 
agreements between the MCSM and FCI results. The detailed comparisons among 
the MCSM, FCI and NCFC methods can be found in Ref. HSj). 

By exploiting the recent development in the computation of the Hamiltonian 
matrix elements between non-orthogonal Slater determinants^^' and the technique 
of energy- variance extrapolation,^ the observables give good agreement between the 
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MCSM and FCI results in the p-shell nuclei. From the benchmark comparison, the 
no-core MCSM is now verified in the application to the ab initio no-core calculations 
for light nuclei. Moreover the application of the no-core MCSM to heavier nuclei is 
expected in the near future. 

3.3. Analysis of intrinsic state 

While ab initio approaches have been studied intensively in light nuclei, it is 
relatively difficult to study the cluster structure in an ab initio way. Among these 
approaches, the Green's Function Monte Carlo first provided the two-a structure 
of the ^Be ground state illustratively.!^ This study has shown the possibility that 
the cluster structure can appear in ^Be, without assuming any cluster structure in 
advance. Generalizing this result, it may be possible to treat cluster structure from a 
pure single particle picture. In this subsection, we show how to visualize the cluster 
state in the no-core MCSM calculation and by analyzing the calculations we discuss 
the appearance of a cluster structure. It is also suited to clarify the relation between 
the shell-model and cluster pictured^ from the shell-model point of view. This view 
point has not been investigated very well yet. Recently, the density profile in the 
lithium isotopes has been investigated by the NCFClSSJ The method has shown how 
to calculate the translationally-invariant density. In Li isotopes, the shape distortion 
and cluster-like structure has been found. Thus, the study of cluster structure has 
become a realistic subject by using the shell-model calculation. 

To extract the cluster structure from the no-core MCSM, we define the intrinsic 
state to visualize the cluster shape in the intrinsic framework which is extracted 
from the angular- momentum-projected wave function. The wave function of the no- 
core MCSM, which is defined in Eq. ()2-5p . is represented as an angular- momentum 
projection of a linear combination of basis states such as 

|^) = P^|^) \^) = Y^fM, (3-1) 

n 

where the total / is assumed to be zero and i^-quantum number and parity projec- 
tions are omitted for simplicity. This linear combination of the unprojected basis 
states, 1^), cannot be considered as an intrinsic state because the principal axis of a 
basis state, is not in the same direction as that of another basis state. There- 
fore we rotate each basis state so that it has a diagonalized quadrupole- moment; 
Qzz > Qyy > Qxx and Qij = 0, (i 7^ j), respectively, following the concept of Ref. 

As a result, these rotated basis states have a large overlap with each other and 
make a distinct principal axis toward the z-axis. The intrinsic wave function |<^™*^) 
is defined as 

^ = (3.2) 

n n 

where the R{{ln) is a rotation operator with Euler's angle Qn- The f2n is deter- 
mined so that the transformed basis state \4>^) = R{f2n)\4'n) has the diagonalized 
quadrupole-moment. The transformed coefficient (by R{Qn)) is derived by the 
relation in Ref. I12p . This state exactly has the same energy after the angular mo- 
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Fig. 8. *Be proton density for |#) (left panels) and intrinsic (right panels) states for various A^;, 
and sliced along the yz plane. The number of basis states is A^'t = 1, 10 and 100 for the lower, 
middle and upper figures, respectively. The slice along the yz plane is the a; = fm plane (left) 
or a; = 1 fm plane (right) for each panel. The size of each box is 8 fm x 8 fm. 



mentum projection. We calculate the one-body density of the intrinsic state such 
as 

pintr(^) = (<?i°f| ^5(r - ri)\^''^''), (3-3) 

i 

where rj denotes the position of the i-th nucleon. 

As an illustrative example, we show the ^Be density in Nsheii = 4 and Hu = 20 
MeV with the JISP16 interaction for J = 0^ states. The Coulomb interaction and 
the contamination of spurious center-of-mass motion are neglected for simplicity. We 
show the proton density (a half of the total density) of the \^) and the intrinsic-state 
density, p™*'', in Fig. El 

The number of basis states is A'^;, = 1, 10 and 100 for the lower, middle and 
upper rows, respectively. The energy is almost converged at Nf, = 100. Each density 
distribution is shown along the yz planes at x = fm and at x = 1 fm. As shown in 
the N), = 1 results, clear deformation and the neck structure to be called a dumbbell 
shape appear. We can see that as the number of basis states, Nh, increases the 
density of |^) are much vague and becomes ordinary prolate rather than dumbbell- 
like because of the mixture of different directions of principal axes of the basis states. 
On the other hand, the intrinsic density has clearer dumbbell-like structure for each 
Nh. In addition, the density distribution of the intrinsic state is almost unchanged 
with respect to Nb. This result indicates the appearance of cluster structure in the 
no-core MCSM. We also check how the cluster shape differs between Ng^f^u = 3 and 
^sheii = 4. We find that the neck of dumbbell shape is more enhanced in Nsheii = 4 
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than in Nsheii = 3. Since the weights of distribution for both sides of the principal 
axis are almost the same, this cluster can be considered as two a clusters. The 
stability of the a cluster is confirmed with respect to Nf, and A^^/ieZZ- This result is 
consistent with the result of the Green's Function Monte Carlo!® With the use of 
this method to draw the density, we can study the appearance of cluster structure 
directly not only for N = Z nuclei but also for the neutron-rich nuclei in the ab 
initio approach. The study of exotic structure including unstable nuclei in the p- 
shell region is in progress. 

§4. Application to neutron-rich Cr and Ni isotopes 

In this section, we discuss the application of the MCSM to the large-scale shell- 
model calculations about neutron-rich Cr and Ni isotopes as examples. We take a 
model space as the pfgQd^ shell, which consists of the Oflp shell, the Ogg/2 orbit, and 
the 1^5/2 orbit. By using such a sufficiently large model space, we aim at a unified 
description of medium- heavy nuclei and at studying the shell evolutioiPSJ'CSli'EBiET) 
and the magicity of = 28, 40, 50 microscopically. 

4.1. Ni isotopes and magicity of N = 28,40,50. 

The nuclear shell structure evolves in neutron-rich nuclei and the magic numbers 
of unstable nuclei are different from those of stable nuclei. The large excitation energy 
of the 2+ yrast state and the small B{E2; 0+ 2+) value in ^^Ni {Z = 28, N = 40) 
might indicate that ^^Ni is a double-magic nucleus, although = 40 is a magic 
number of the harmonic oscillator, not a magic number of the nuclear shell model. 
On the other hand, the small excitation energies of the 2^ yrast state and the large 
B{E2;0'^ — )• 2"*") values in Cr {Z = 24) isotopes of ~ 40 suggest rather strong 
deformation. This change of the = 40 gap has been studied theoretically.'^ ^^Ni, 
which has Z = 28 and A^ = 50 doubly magic numbers, has also been investigated to 
discuss its magicity and the size of A^ = 50 gap^ 

In the sd-shell and the light p/-shell regions, we can describe properties of stable 
nuclei in relatively small model spaces. However, we sometimes need a large model 
space to describe the properties of neutron-rich nuclei. In order to discuss neutron- 
rich Ni isotopes up to A^ = 50, it is essential to include the effects of excitation 
across the Z = 28 and A^ = 50 gaps by adopting the pfggd^ model space. Concern- 
ing this model space, M. Honma et al. proposed the A3DA effective interactiorPSJ 
which consists of the GXPFIA,®^' JUN45,'^ and G-matrix effective interactions 
with phenomenological modifications. It has succeeded in describing the neutron- 
rich Cr and Ni isotopes under a severe truncation of the model space utilizing the 
few-dimensional basis approximation.'^ In this work, we use the new version of the 
MCSM method, which enables us to precisely evaluate the exact shell-model energy 
without any truncation and discuss the effective interaction. 

4.2. Effective interaction for pfggd^ shell 

In this section, we discuss the A3DA effective interactioiJ® and its improve- 
ment. The two-body matrix elements (TBMEs) of the A3DA interaction consist of 



20 



three parts. The TBMEs of the pf shell are those of the GXPFIA interaction]^ 
which is successful for describing spectroscopic properties of light p/-shell nuclei. 
The TBMEs of the f^pgg shell related to the 0^9/2 orbit are those of the JUN45 
interaction.lSJ The GXPFIA and JUN45 interactions were determined by combin- 
ing microscopically derived interactions (G matrix) with a minor empirical fit so as 
to reproduce experimental data. The other TBMEs are from the G-matrix effec- 
tive interaction}^ which is calculated from the chiral N3L0 interaction.HSJ The 
Coulomb interaction is not included and the isospin symmetry is conserved. The 
G matrix is calculated for the pfsdg shell with ^"^Ca as an inert core and the core- 
polarization correction is included perturbatively. The single-particle energies and 
the monopole interaction are adjusted to reproduce the GXPFIA and JUN45 predic- 
tions for the pf shell and gg/2 orbits, and the Woods-Saxon single-particle energies 
of stable semi-magic nuclei for the other part. 

The original A3DA interaction failed to describe some nuclei around N ~ 40. 
We modify mainly single-particle energies and monopole components related to the 
0(79/2 orbit by comparing the results of the calculations with the experiments. These 
calculations are far beyond the current limitation of the conventional diagonalization 
method, and the MCSM method enables us to perform this comparison. 

4.3. MCSM results of the neutron-rich Cr and Ni isotopes 

We performed systematic calculations of the 0"*" and yrast states of neutron- 
rich Cr and Ni even-even isotopes using the MCSM method and the modified A3DA 
interaction. We took 50 basis states for the MCSM with the refinement procedure, 
which is discussed in Sect l2.2l The energies were extrapolated by the energy-variance 
extrapolation method and the other values were not. The effective charges are taken 
as (ep,e„) = (1.5,0.5)e. 
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Fig. 9. The excitation energies of 2]^ states (left) and B(E2;0f — >■ 2f) values (right) obtained by 
the MCSM for Cr isotopes. Experimental data are taken from Refs.'^''^ 

Figure shows the 2"*" excitation energies and the B{E2) transition probabilities 
of neutron-rich Cr isotopes. The MCSM results well reproduce the experimental 
values while the modest overestimation remains. The Cr isotopes do not show any 
feature of = 40 magicity, while the characteristics of A^ = 28 magicity can be seen. 
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namely, a sudden increase of excitation energy and slight decrease of the B{E2) 
value. On the neutron-rich side, the excitation energy decreases and the B[E2) 
value increases gradually as the neutron number increases, which implies gradual 
enhancement of the quadrupole deformation. 

Figure [TOl shows 2^ excitation energies of Ni {Z = 28) even-even isotopes from 
^^Ni to '^'^Ni. The large 2+ excitation energy of ^"^Ni {N = 28) indicates Z = 28, 
= 28 double magicity. The large value of the calculated 2"'" excitation energy of 
'^^Ni (A^ = 50) suggests Z = 28, = 50 double magicity. The large 2^ excitation 
energy of ^^Ni (A^ = 40) indicates A^ = 40 magicity. The calculated values reproduce 
the experimental values well. 

Figure [TT] shows B{E2;0'^ — ?> 2^) for neutron-rich Ni isotopes. The small value 
of B{E2; 0^ — )• 2+) at A^ = 40 indicates A^ = 40 magicity. Neither the theoretical nor 
experimental value of B{E2; 0+ 2+) at A^ = 28 is smah unhke that at the A^ = 40, 
and the theoretical B{E2;0'^ 2+) value at A^ = 50 becomes large in comparison 
with those of neighboring nuclei. It suggests that at A^ = 28, 50 magicity is broken to 
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Fig. 14. Total energy surface of 0^ (left) and (right) states of "^Ni. The positions of red circles 
represent quadrupole deformations of the MCSM basis states before projection. The areas 
of those circles represent the overlap probabilities of the basis states and the resulting wave 
function. 

some extent for ^^'^^Ni, respectively. Figure [T^ shows the occupation number of the 
neutron (79/2 orbit. The occupation numbers of 0^ and 2^ states are very close for 
Ni isotopes besides ^^''^^Ni {N = 40, 50). The occupation numbers of the states of 
68'78]\jj show a breakdown of the closed-shell structure. Figure [13] shows two- neutron 
separation energies. The calculated values of neutron-rich nuclei are smaller than 
experimental values. This means that the binding energies of neutron-rich nuclei 
are underestimated. The values of S2n increase slightly by considering the Coulomb 
energy, but calculated values are still smaller than experimental values. 

Figure [T3] shows the total energy surface of ^^Ni provided by the Q-constrained 
Hartree-Fock calculation.'^ There are three minimum points for ^^Ni. Figure [HI also 
shows quadrupole deformations of the MCSM wave functions of the Of and 0^ states. 
The scattered circles correspond to the basis states in the MCSM wave function. 
The position of the circle indicates the quadrupole deformation of the basis state 
before projection. The area of the circle is proportional to the overlap probability 
of the projected basis and the resulting wave function. It is quite clear that the 0^ 
state of ^^Ni corresponds to a spherical shape and the 0^ state corresponds to an 
oblate shape. Spherical and oblate components are mixed to some extent, but the 
components of the prolate minimum hardly mix. 

Furthermore, we consider the magicity and the energy gaps for Ni isotopes by 
using the effective single particle energies (ESPEs).E3 Fi gure 1151 shows ESPEs of 
the neutron orbits for Ni isotopes. The /7/2-P3/2 S^-P at = 28 is 7.1 MeV and 
gives the magicity to ^^Ni. The ^9/2-^5/2 g^-P at = 50 is 4.2 MeV and gives the 
magicity to ^^Ni. This is partly due to the additional lowering of the 39/2 orbit 
caused by pairing correlation between two neutrons in the gg/2 orbit, and also due 
to the effect of the two-neutron repulsive monopole interaction originating in the 
three-nucleon force like in exotic oxygen isotopes .'22' The ^1/2-59/2 g^'P at = 40 is 
2.6 MeV, which is smaller than the N = 28, 50 gaps. Figure [TBI shows the ESPEs of 
the neutron orbits for = 40 isotones. As the proton number of increases from 
Z = 20 to Z = 28, the ESPE of lowers and the A^ = 40 gap becomes larger. 
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Fig. 15. ESPEs of the neutron orbits for Ni 
{Z — 28) isotopes. 
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Fig. 16. ESPEs of the neutron orbits for A'' = 
40 isotones. 



Because of this evolution of the = 40 gap, the properties of ~ 40 nuclei depend 
on the proton number. 




Fig. 17. ESPEs of the proton orbits for Ni isotopes (left) and for A'^ = 40 isotones (right). 



In Fig. [T71 the ESPEs of the proton orbits for Ni isotopes and for = 40 
isotones are shown. In the former, rapid lowering of the orbit from = 40 to 
50 is clearly seen as suggested in,'i3''23 while narrowing of Z = 28 gap is also visible 
there. Such changes are responsible partly for the origins of the structure evolution 
in these Ni isotopes. 

§5. Summary and future perspectives 

We have developed a new generation of the MCSM by introducing the conjugate 
gradient method and the energy-variance extrapolation, which enhance the applica- 
bility of the MCSM greatly. We have two major scopes of this framework: ab initio 
shell-model calculations and conventional shell-model calculations assuming an in- 
ert core. In the former, we have compared the MCSM results with the exact FCI 
calculations to demonstrate the validity of the MCSM framework and its feasibil- 
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ity beyond the limit of the FCI in Secti3l In addition, we have proposed a novel 
method to discuss the intrinsic structure and demonstrated that the cluster struc- 
ture appears in shell-model-type calculations based on the harmonic-oscillator-basis 
wave function. In the latter, we discussed in SectSlthat the MCSM enables us to 
perform shell-model calculations of neutron-rich Cr and Ni isotopes in the pfggd^ 
model space in which the isospin symmetry is conserved. We proposed a "modified 
A3DA" interaction which reproduces the low-lying spectra of neutron-rich Cr and Ni 
isotopes and guides us towards a unified description including ^^Ni, ^^Ni and ''^Ni, 
with magic numbers 28, 40, 50, respectively. The prediction of ^^Ni is especially 
interesting to see the evolution of shell structure. On the other hand, Cr isotopes 
do not show any feature of = 40 magicity and the collectivity enhances as the 
neutron number increases. The MCSM and newly proposed effective interaction are 
expected to provide us with a unified description of p/-shell nuclei. 

The current status of the computer-code development was also reported in 
Sect 1231 At the present stage, we have obtained good parallel scalability of our 
code up to 10^ CPU cores via early access to the K computer at RIKEN AICS^^' as 
measured by the benchmark test. However, such good scalability is not always ob- 
tained and further development is in progress. This activity is promoted strongly as 
a part of the activities of HPCI Strategic Programs for Innovative Research (SPIRE) 
Field 5 "The origin of matter and the universe". 

By utilizing both the developed code and the K computer, we promote further 
large-scale shell-model calculations as a part of the SPIRE activities. We plan to 
perform systematic study with ab initio calculations of light nuclei in Nsheii = 5 
and some states in Nsheii = 6. Concerning the medium-heavy nuclei, because it is 
difficult to cover whole region of the nuclear chart, we will choose some interesting 
nuclides as subjects of our investigation, and will perform shell-model calculations of 
these nuclides with the two-major-shell model space. For example, the shell-model 
calculations of ^'^'^Te, ^■^'^Te, and ^^'^Nd are extremely interesting to study double 
beta decay and the nuclear matrix element of neutrinoless decay. We also continue 
to study the systematic calculations of neutron-rich p/-shell nuclei to discuss the 
shell-evolution phenomenon. 
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Appendix A 

Numeration with projected Slater determinants 

In this appendix, we show some equations which are needed to perform the 
calculation discussed in SectlJJ 

At the beginning, we define a deformed Slater determinant, 

k=i \i=i J 

which is parametrized by the complex A'^gp x A^f matrix D with the normalization 
condition D^D = 1. Nf and Asp are the numbers of fermions and single-particle 
states, respectively. The |— ) denotes an inert core in the conventional shell-model 
calculations or the vacuum in ab initio shell-model calculations. Because we do 
not mix the proton and neutron space in practical calculations, the wave function 
is written as a product of proton and neutron Slater determinants, namely, \4>) = 
|<^proton) "Si | (/'neutron ) • For simplicity, we do not write this isospin degree of freedom 
explicitly. One can easily reproduce the equations representing the explicit proton- 
neutron degree of freedom by taking D of the proton-neutron sector as zero such 
as 

■ D"" 



where and represent Slater determinants of protons and neutrons, respec- 
tively. 

The angular- momentum, parity projector Pj^j^ in Eq. ()2-6p is performed by dis- 
cretizing the integral concerning the Euler angles such as 

A 

where the A denotes an index of mesh point of the discretization (here, a set of Euler's 
angle Q = (a, /3, 7) and parity variable vr'** = ±1). In this paper, the parity projection 
is described by the summation of 2 mesh points such as P"^ = ^-^^ = X]a=i t^^^^ n^^^ 
with vr^^) = ^, TT^^^ = ^, n^^^ = 1, and U^"^^ = U with 77 being the parity-conversion 
operator. W^'^^ is a weight of the mesh point A, and R^^^ is a product of the rotation 
and parity-conversion operators such as 

Kfic'^ = ^I?ik(«A,/?A,7A)^('\ 

^(A) _ ^iaxJz^il3\Jy^i'y\Jz]j{\) _ (A-4) 

Note that the operator 7?^ does not change the form of a Slater determinant, i.e., 

|<AW) = 7?W|<A„), (A-5) 

where a matrix 7)"(^) represents the single Slater determinant \(f)^^), thanks to the 
Baker-Hausdorff 's theorem.^^^ 
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The norm matrix and hamiltonian matrix spanned by N Slater determinants 
are written as 

J^mM,nK = {4'm\PMK\4>n) (A-6) 
T~LmM,nK = {(t>m \HPilM. (A-7) 



The coefficient, fnK in Eq. ()2-5p . is determined by solving the generalized eigenvalue 
problem 

y^^T-l-mM,nKfnK = £ y^^J^mM,nK fnK , (A'S) 
nK riK 

and the normalization condition (if']??') = 1. The lowest eigenvalue of £ is taken as 
En if you would like to obtain the yrast state. 

By combining Eqs. ()A-3P and ()A-6p . the norm matrix is calculated as 

A A 

with 

(0^|4^)) = det(z)™ti5"W). (A-10) 
In the same way, the hamiltonian matrix is obtained as 

nmM,nK = Y^^MK^^rnlH R'^^^ {cPn) (ATI) 
A 

where the generalized density matrix, p^^\ and the self-consistent field, are 
defined as 

^(A) ^ ('^"-14^-1^"' ) ^ pn(A)pMz)n(A))-l^M).. (A.12) 



4'^=Ewl-^ (A-13) 

jl 

with Vijki = Vijki — Vijik- The trivial summations and their indices for the matrix 
products are omitted for readability. The indices m, n of p(^) and r(^) are also 
omitted. 

The most-time-consuming part is the calculation of the r-^\ which can be 
rewritten following the idea of Ref. 

rj^'^ = Y.^a,pi'^ (AT4) 

b 

where a = {i,k), and b = Because Vab is a block-antidiagonal form owing to 

the symmetry of the Hamiltonian, Eq. ()AT4p is calculated as the products of the 
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dense block matrices and the dense matrices in terms of the indices a, b, A efficiently 
avoiding trivial zero matrix elements of Vab- This method is referred to as the matrix- 
matrix method in RefJ42p. This matrix-matrix method enables us to use a CPU 
utilizing the BLAS level 3 library quite efficiently, and the performance reaches 
70 ~ 80% of the theoretical peak performance.^^' 

This efficient computation of F^^^ is useful also for the evaluation of the energy 
gradient, which is essential for the conjugate gradient method. The energy gradient 
of the Slater-determinant coefficients is written as 

^ = (l-D^D^^) /4M/ni.<l'^0m|4'^) (A-15) 

M,n,K,X 

1 - + rW) + (Tr((t + ^rW)pW) - i?^)) pWz)-. 

To evaluate the energy variance {AH'^)n = {H'^)n — Ej^, the expectation value 
of the H"^ with the wave function in Eq. ()2-5p is written as 

mH'\^N)= E fUifnKWi;j,^\cl>^\H'\<l>i^^). (A.16) 

m,M,n,K,X 

From Ref|17l) , the matrix element of the Hamiltonian squared is computed such as 



(A)^ - '"ijki&kiLp'^c.fi^sA\s\- (A- 17) 

i<j,k<l,a<l3,j<5 



2 



+Tr((f + rW)(i - pW)(t + rW)pW) + (^TY(pW(t + irW)) 

/I (A) _ ^(A) JA) (A) (A) . 
A-ijkl-PikPjl - Pil Pjk (A-18j 

= (1 - P^'^Ml - P'-^I - (1 - p'-^.lil - P^'^hk. (A-19) 

The most-time-consuming part in the evaluation of the energy variance is to calculate 
the first term of the right-hand side of Eq. ()A-17p . By substituting {i,j), {k, I), (a, /3), 
and (7, 5) by a, b, c, and d, respectively, this term is efficiently calculated as the 
products of the matrices, namely, X]a6cd '^ab&bcVcdAda- Note that the Vab has a block- 
diagonal form, which again enables us to use the BLAS level 3 library, avoiding trivial 
zero matrix elements. 

Another formulation to compute the expectation values in projected Slater deter- 
minants can be found in Refs.,'2Sl''EHli in which the two-body interaction is decomposed 
into a sum of the squares of the one-body operators. 
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